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Abstract. The three-dimensional, hydrodynamic stability of the solar tachocline is investigated based on a rotation 
profile as a function of both latitude and radius. By varying the amplitude of the latitudinal differential rotation, we 
find linear stability limits at various Reynolds numbers by numerical computations. We repeated the computations 
with different latitudinal and radial dependences of the angular velocity. The stability limits are all higher than 
those previously found from two-dimensional approximations and higher than the shear expected in the Sun. It is 
concluded that any part of the tachocline which is radiative is hydrodynamically stable against small perturbations. 



1. Motivation 

Hclioseismology has provided us with information on the 
internal rotation of the Sun. The angular velocity depends 
on latitude as well as radius. The dependence is mainly a 
latitudinal in the bulk of the convection zone, whereas the 
solar radiative core rotates nearly uniform with an angular 
velocity of the convection zone at about 30° latitude. The 
transition from the differential rotation in the convection 
zone to the uniform rotation in the core is thin - probably 
thinner than 5% of the solar radius - and is called the 
tachocline. 

The convection zone is thermally overcritical and the 
stability of shear flows is not an issue there. The shear in 
the tachocline however is latitudinal and radial and may 
be subject to shear-flow instabilities. If the tachocline is 
turned into a turbulent layer, a problem arises with the 
mixing of elements, most notably with Lithium which will 
be destroyed in nuclear fusion up to 0.68 solar radii (i?©), 
just below the convection zone. A turbulent tachocline 
mixing the Lithium into its fusion zone would contradict 
the observed relatively high Lithium abundance at the 
surface of the Sun. 

The situation in the solar tachocline was described by 
Spiegel & Zahn (1992) as exhibiting horizontal turbulence, 
excited by hydrodynamic shear-flow instability. It was ar- 
gued that the two-dimensional turbulence provides angu- 
lar momentum transport but inhibits too strong mixing 
of Lithium below 0.68i?©. 

The hydrodynamic stability of the latitudinal shear in 
the tachocline was studied by Watson (1981). The depen- 
dence on colatitude 9 was fi = fi eq (l — 012 cos 2 8), where 
Sl e q is the equatorial angular velocity at the bottom of 
the convection zone and a-i is a parameter for the rel- 
ative equator-to-pole difference of the angular velocity. 



Watson found instability for a differential rotation with 
a2 > 0.286. At that time, this was supposed to be in the 
vicinity of the value of solar differential rotation, and the 
result was not definitely deciding between a turbulent and 
a stable tachocline. The investigation was two-dimensional 
arguing that the stable stratification will not allow signifi- 
cant radial flows, but thereby the radial shear is neglected, 
too. 

Charbonneau et al. (1999a) extended the linear stabil- 
ity analysis to various rotation profiles of the form 

fi = n eq (l - a 2 cos 2 6> - a 4 cos 4 6>). (1) 

Different layers of the tachocline were associated with dif- 
ferent results from helioseismology in terms of a 2 and 
Q!4. The upper layer which is thought to be penetrated 
by convective overshooting was found to be unstable to 
the shear, whereas the lower layer with smaller latitudinal 
shear turned out to be stable. In a step towards the three- 
dimensional stability of the tachocline, Dikpati & Oilman 
(2001) studied the two-dimensional system allowing for 
deformations in the third - the radial - dimension. As a 
result, the critical differential rotation for instability was 
reduced to 11% in the overshoot part of the tachocline. 

In order to address the entire tachocline and since the 
tachocline is a place where latitudinal and radial shear 
meet, we investigate the stability of the three-dimensional 
rotation profile. Although it is very reasonable to assume 
that radial flows will be weak, the variation of the latitu- 
dinal differential rotation with radius across weakly cou- 
pled spherical layers could provide different results for the 
stability of the tachocline. We will briefly summarize the 
numerical background in Section 2, provide details of the 
computational results in Section 3 with different rotata- 
tion profiles, and summarize our findings in Section 4. 
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Fig. 1. Vertical cross-section through the solar tachocline 
with contours of the assumed angular velocity depending 
on radius and latitude as given by (2). Highest angular 
velocity is at the equator (bottom right), and radially con- 
stant rotation occurs at a hcliographic latitude of 30°. 

We find that the tachocline is hydrodynamically stable in 
all the configurations studied. The paper is thus intended 
to be a preparatory step towards fully three-dimensional 
studies of the MHD stability of the tachocline. 

2. Computational setup 

The rotational profile depends on both latitude and radius 
in this study. Between the inner and outer radius of the 
tachocline, n — 0.65 and r = 0.7 respectively, the angular 
velocity is defined by 



n(r, 9) = a 



cq 





\r a -r 









,(2) 



where 9 denotes the colatitude in the spherical shell, r is 
the radial coordinate, and f2 eq is the equatorial angular 
velocity. The profile implies that the rotation velocity at 
the inner boundary of the computational domain, r = fj, 
is the one at 9 = 60° (or 30° heliographic latitude). This 
appears to be a valid assumption in agreement with var- 
ious helioseismological inversions (recently e.g. Schou et 
al. 2002), whereas core rotation may be adopted at larger 
latitudes higher up in the convection zone. The formation 
of the tachocline rotation profile is supposed to be caused 
by rotating convection on top of it and reduction of dif- 
ferential rotation by magnetic fields in the solar core at 
the bottom. Such a profile causes a meridional circulation 
reaching a steady-state in competition with the aforemen- 
tioned (or other tachocline-forming) effects (see e.g. Sule 
et al. 2005). In our linear analysis, this flow has no effect on 
the stability of the non-axisymmetric modes investigated 
in this Paper. 

We employ the incompressible, viscous Navier-Stokes 
equation and linearize the problem. We can then separate 



the axisymmetric background rotation U from the non- 
axisymmetric flow u. The latter is evolved by numerical 
computations. The normalised equation of motion reads 



du 

— — = uxVxU+UxVxu 
at 

-Vp- V(« ■ U) + Am, 



(3) 



and the continuity equation V • u = holds. The equation 
is evolved with the spectral spherical code by Hollerbach 
(2000). We look for exponentially decaying or growing so- 
lutions in order to find the critical differential rotation 
a 2 of the marginal case. The actual integration employs 
the radial components of the curl and the curl-curl of the 
equation, thereby eliminating the gradient terms. The lin- 
ear (viscous) part of the full equation of motion is evolved 
implicitly, while the nonlinear parts are integrated explic- 
itly with the advection term and forces being computed 
in real space. We kept this splitting even in our linearized 
problem, since a fully implicit scheme would have required 
a substantial modification of the code. The Am is thus 
treated implicitly, the other rhs terms (two remaining af- 
ter curling) are computed in real space and are used for a 
second-order Runge-Kutta integration. 

The normalisation of the equation with the viscous 
time r u = r\jv and the length scale r leads to the 
Reynolds number 



Re 



1 o"cq 



(4) 



as a free parameter which is essentially a variation of the 
viscosity v since radius and f2 Gq are sufficiently well known. 
The solar Reynolds number in the tachocline is - in terms 
of the definition of Q - about 10 14 . We try to achieve 
time series for numerically demanding Re > 10 4 . By com- 
parison with known results from inviscid two-dimensional 
analyses, we find that the critical viscous differential ro- 
tation at Re > 10 3 or 10 4 is already sufficiently close to 
the inviscid value. 

The velocity is decomposed into toroidal and poloidal 
potentials, u = V x (er) + V x V x (fr), where f is the 
unit vector in radial direction. The potentials are again de- 
composed into radial Chebyshev polynomials and spher- 
ical harmonics. Since it is the potentials being evolved, 
the continuity equation is fulfilled automatically. The den- 
sity is constant throughout the computational domain. In 
a thin shell of 5% of the solar radius, this is a reason- 
ably good approximation of the true situation. We also do 
not take any deformation of the tachocline into account. 
Top and bottom radius of the tachocline are constant 
over latitude. Charbonneau et al. (1999b) found a pro- 
late tachocline whose equatorial part is located at slightly 
smaller radius than the polar end. We assume that the 
difference of 3.5% in location of the tachocline has negli- 
gible effect on the results, but may be an issue of future 
investigations. 

The radial boundary conditions for the velocity per- 
turbations are stress- free at both r\ and r Q . At Reynolds 
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Fig. 2. Streamlines of the symmetric eigen function of the 
computation with a three-dimensional profile Q(r,8) on 
the surface of the tachocline at r = 0.7. Two circulation 
cells are found on each hemisphere. 

numbers of Re > 10 4 , high spectral resolution was neces- 
sary to obtain reliable results. Up to 80 Chebyshev and 80 
Legendre polynomials were used to resolve the flow prop- 
erly. 

The azimuthal modes of the problem described by Q 
are decoupled, and we can study the stability of individual 
m-modes separately. Moreover, even and odd latitudinal 
modes (symmetric and antisymmetric modes with respect 
to the equator) decouple, and we will have a look into the 
critical differential rotation for the excitation of instability 
of the two kinds separately. The radial modes all couple 
and do not provide results on the stability of individual 
radial wavelengths. 

3. Results 

3.1. Stability of various solutions 

In a set of fiducial computations, we applied a purely lat- 
itudinal profile of the angular velocity. An m = 1 mode 
is evolved with the profile of (1) where = and a.i 
is varied. Since we solve a viscous problem, the critical 
differential rotation, a™', depends on the Reynolds num- 
ber. The result is already very close to Watson's inviscid 
result for Re > 1000. This is good reason to assume that 
the numerical solutions reaching Re ~ 30 000 are suitable 
approximations for the solar plasma. 

Despite allowing for radial motions, the evolution pro- 
vides solutions which are nearly toroidal and do not show 
significant radial flows. They are surface flows forming 
two cells on each hemisphere with stream lines through 
the poles. Figure shows a representation of the flow in 
the spherical surface. The graph has to assume that the 
poloidal component of the velocity is zero though, which 
is not entirely true. 

The second step involved the rotation profile of (2) 
for which the critical steepness of the differential rota- 
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Fig. 3. Lines of marginal stability for the combined lati- 
tudinal and radial shear (solid line) and the purely latitu- 
dinal shear (dashed line). Differential rotation denotes the 
percentage by which the pole's angular velocity is slower 
compared with the equatorial one. It is expressed by a,i 
from (1) and (2) here. 
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Fig. 4. Lines of marginal stability for the antisymmetric 
m = 1 solutions caused by the combined latitudinal and 
radial shear (solid line) and the purely latitudinal shear 
(dashed line). 

tion, a™*, is again sought for various Reynolds numbers. 
Figure 13 shows the lines of marginal stability, i.e. the crit- 
ical differential rotation, versus Reynolds number for the 
symmetric m — 1 mode. The solid line refers to profile (2), 
the dashed line is the latitudinal profile and converges to 
the result by Watson (1981) for Re — •> oo. 

The most easily excited patterns of m — 1 are always 
symmetric with respect to the equator. We can also look 
for the stability of antisymmetric patterns and find the 
results shown in Fig. They are more stable than the 
symmetric configurations with an f2(r, 9) profile. The anti- 
symmetric solutions from the Q(9) profile have also higher 
critical differential rotation values than their symmetric 
counterparts. 

The patterns drift with a certain velocity in azimuthal 
direction around the solar axis. Since the equations hold 
for the nonrotating system, we can directly convert the 
pattern rotation into physical times. The pattern rotation 
period for the two-dimensional and the three-dimensional 
profile of f2 is shown in Fig. [5] by a dashed and a solid 
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Fig. 5. Rotation period of the flow pattern for the 
two-dimensional (long-dashed) and the three-dimensional 
(solid) rotation profile. The numbers are computed as- 
suming an equatorial rotation period of 25.44 d (fi q = 
455 nHz) which is plotted with a dash-dot line. Periods 
are computed at marginal stability; the polar rotation 
period for this critical differential rotation is plotted as 
short-dashed and dotted lines for the 2D and 3D cases, 
respectively. 



line, respectively. The actual solar rotation periods are 
also plotted. The reference equatorial period of 25.44 d 
(fi e q = 455 nHz) is given as dash-dot line. The pattern ro- 
tation periods are determined at marginal stability. Since 
the marginal case gives us a value for the differential ro- 
tation, we also plot the polar rotation period P po \ c versus 
Re. The short-dash line is for the two-dimensional 0,(9) 
profile, the dotted line is for the three-dimensional case 
described by (2). 

The pattern rotation periods are always between the 
equatorial and polar rotation periods, in agreement with 
the 2D results by Charbonneau et al. (1999a). While the 
patterns from the 2D-S1 profile are close to the polar ro- 
tation period, the patterns for the 3D profile rotate with 
nearly the average rotation period between the polar and 
equatorial ones. 

We can compute the time after which the pattern is 
passed by a given point on the equator. This time is often 
called lap time. Assuming an equatorial rotation period of 
25.44 d (fi e q = 455 nHz), we find a lap time of 91 d for 
the 2D case, and a lap time of 78 d for the 3D case. 

Modes with higher azimuthal mode numbers require 
significantly higher differential rotation for instability. The 
antisymmetric m = 2 mode which is symmetric with re- 
spect to the equator was found to be stable even in the 
entire parameter range covered by Fig.|3 This is in agree- 
ment with the inviscid, two-dimensional stability analysis 
by Charbonneau et al. (1999a). The stability lines for the 
antisymmetric m = 2 mode are shown in Fig. El We could 
not find instability for any m — 3 mode in the range cov- 
ered by Fig. 
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Fig. 6. Lines of marginal stability for the antisymmetric 
m = 2 mode resulting from the combined latitudinal and 
radial shear (solid line) and the purely latitudinal shear 
(dashed line). 



3.2. Effects of buoyancy 

Little influence is expected from the stable temperature 
gradient in the tachocline. Since we do have the chance 
to prove this in our three-dimensional simulations, we 
demonstrate the effect of a negative buoyancy force on 
the stability of the differential rotation. The Navier-Stokes 
equation is extended by the buoyancy force and reads in 
the non-dimensional form: 



du 
~dt 

dQ 
~dt 



uxVxU+UxWxu 
+Ra6r - Vp - V(tt ■ U) + Am, 

U VB-u- VT+— A6, 
Pr 



with a background temperature profile of 



T : 



n 



(5) 
(6) 

(7) 



and the Prandtl number being the ratio between viscosity 
and thermal diffusivity, Pr = vj\- The Rayleigh number 
in © is 



Ra = 



9 a(Ti - T Q )rl 



(8) 



where g is the gravitational acceleration, a is the coeffi- 
cient of volume expansion, and Tj and T Q are the tempera- 
tures at the two boundaries. In the Boussinesq formulation 
used here, the presence of a sub-adiabatic temperature 
gradient actually translates into a negative value of Ra. 
We set our version of the Rayleigh number to a value as 
small ("as negative") as Ra = —10 s in order to see any 
notable effect on the flow. The Prandtl number is set to 
unity. 

The critical differential rotation for a growing symmet- 
ric m = 1 mode at Re = 10 4 increases slightly to 53.3% 
as compared with the non-buoyant value of 52%. This is 
in line with the fact that the solutions contain nearly hor- 
izontal motions. 
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Fig. 7. Contours of the assumed angular velocity with a 
radial dependence as in (|1(J|) . Highest angular velocity is 
again at the equator (bottom right), and radially constant 
rotation occurs at a heliographic latitude of 30°. 

3.3. Effects of higher-degree terms 

The differential rotation has been expressed by a cos 2 9 
dependence in latitude and a linear r dependence over 
radius. We are studying the stability of the three- 
dimensional tachocline setup with higher-degree depen- 
dences in this Section, such as the cos 4 9 term and a non- 
linear radial dependence of f2. 

Including the cos 4 9 term in (2) yields an angular ve- 
locity profile of the form 

fl(r, 9) = O eq < 1 — a2 cos 2 9 — otA cos 4 9 — ■ 



n 



< i 2 1 — — cos 9 



Ck4 



1 

16 



cos 4 9 



(9) 



The computations for the full profile could not easily be 
extended beyond Re = 5000. But the results for the possi- 
ble Re and for the easier two-dimensional Q(9) show no de- 
creased critical differential rotation when the cos 4 9 term 
is included. This also holds for the extreme case of 
carrying the shear alone (ot2 = 0). 

Starting from @ , a modified Sl-profile was constructed 
in order to find any influence of the particular radial de- 
pendence of the latitudinal shear on the results. We used 



fi(r, 9) 



a, 



1 — Qf2 COS 



| I _ C os 2 0U[ 1 



(10) 



The radial profile introduces two inflection points. It does 
not apply the usual error function for the radial i7-step 
in order to obtain an exact dVt/dr = at both inner 
and outer boundary. A graphic representation is shown in 

Fig.0 
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Fig. 8. Lines of marginal stability for an angular velocity 
profile given by (|10|l having two inflection points over ra- 
dius. The solid line refers to the non-buoyant case, Ra = 0, 
the dashed line has Ra = — 10 8 . 

The critical differential rotation for instability is shown 
versus Reynolds number in Fig. |SJ We were able to reach 
extraordinary Reynolds numbers of 10 5 , probably because 
of the vanishing dil/dr at the boundaries. Instability does 
not occur at reduced differential rotation, and the line of 
marginal stability actually converges to the simpler profile 
(2) for high Re. The additional effect of buoyancy on the 
stability of the profile 11U|) with inflection points is also 
shown as a dashed line. 

4. Summary 

A fully three-dimensional, linear analysis of the stability 
of the solar tachocline was carried out. If radial variation 
of the angular velocity is included in the model, the max- 
imum pole-equator difference of the angular velocity can 
be as large as 52% for a symmetric m = 1 mode, before 
instability sets in. Two-dimensional analyses have deliv- 
ered much lower critical values. The difference of 3D ver- 
sus 2D is not radial flows emerging from the extension in 
the third dimension, but it is changed stability conditions 
emerging from the radial shear and radial dependence of 
the differential rotation. 

Other modes such as higher m or different flow sym- 
metries do not get unstable at lower critical differential 
rotation values under the influence of a three-dimensional 
rotation profile. 

The stabilizing effect of the temperature gradient has 
been added, but since all the unstable modes are very 
nearly horizontal, the influence is small. The assumption 
that horizontal motions dominate is valid even without a 
stabilizing temperature gradient. However the assumption 
that spherical shells of infinitesimal thickness do not in- 
teract with each other is not applicable, according to our 
results. One may argue that the viscosity in the computer 
simulations is much too high, but the variation of the re- 
sults is small at Re > 1000. This is an indication that 
computations with Re = 10 4 or higher are a good approx- 
imation of the near-inviscid solar case. We conclude that 
all parts of the tachocline not being affected by convective 
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overshooting are stable. We believe that the answer on the 
question of the dynamical state of the tachocline will be 
given by MHD computations of a three-dimensional do- 
main. 
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